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Abstract 

Using all atom molecular dynamics simulations, we report spontaneous unzipping and strong binding of 
^ 1 small interfering RNA (siRNA) on graphene. Our dispersion corrected density functional theory based cal- 

culations suggest that nucleosides of RNA have stronger attractive interactions with graphene as compared 
to DNA residues. These stronger interactions force the double stranded siRNA to spontaneously unzip and 
bind to the graphene surface. Unzipping always nucleates at one end of the siRNA and propagates to the 
other end after few base-pairs get unzipped. While both the ends get unzipped, the middle part remains in 
double stranded form because of torsional constraint. Unzipping probability distributions fitted to single 
exponential function give unzipping time (x) of the order of few nanoseconds which decrease exponen- 
tially with temperature. From the temperature variation of unzipping time we estimate the energy barrier to 
unzipping. 
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I. INTRODUCTION 



The interaction of nucleic acids with carbon nanotubes (CNTs) and graphene have attracted 
much attention due to the potential applications in nanotube separation yj, |2D, sensing yi |4j], 



sequencing [5H7J] and nanomedicine |[ Sl4l 
nucleic acids with CNT or graphene 



. Basic understanding of the interaction mechanism of 



1711 is essential for such applications. Small interfering 



ribonucleic acid (siRNA) molecules are a class of double stranded non-coding RNA which are 
typically 21 to 23 nucleotides in length. The properties of siRNA are actively being studied due 
to their potential influence on cell functionality and applications in medicine to achieve RNA 
interference (RNAi) Ill8l - l20n . The mechanism of RNAi involves RNA-induced silencing complex 
(RISC) comprising of Dicer, Argonaute2 and siRNA binding protein that induces unzipping of 
siRNA into two single strand RNAs Il2ll - l25ll . One of these two strands acts as a guiding strand 
to form specific base-pairs with mRNA and silences the gene. For using RNAi technology in 
the treatment of diseases such as cancer, HIV, viral infections and eye diseases efforts are 
being made for the efficient and safe siRNA delivery systems to achieve the desired RNAi effect. 



Dendrimers 1127 



IsJio 



150 are good carriers of siRNA into disease 



infected cell. The discovery of graphene has led to its possible use in efficient delivery of siRNA 
as well as various oligonucleotides. Several attempts have been made to study the properties 
of nucleic acid interaction with graphene 
between graphene and siRNA/DNA has not been understood. 



3 lh . However the interaction 



In this paper, we show the unzipping of siRNA on graphene and subsequent binding between 
the two. Graphene is a free-standing two dimensional monolayer of carbon atoms arranged 
into a honeycomb lattice 11321 . 13311 . Graphene is a partially hydrophobic molecule: its faces are 
highly hydrophobic while edges, depending on functionalization, can be hydrophilic in nature 
113411 . Thus it has the ability to pass through hydrophobic lipid bilayer and can also interact 
with the hydrophilic head groups. Dispersion interaction includin g 71 — 71 stacking serves as 
the major attractive interactions between the non-polar molecules j35\\ . The importance of 
dispersion interactions has been verified recently in analyzing the structure and energetics of the 
graphene-nucleobase complexes [360 and in studying the unzipping of siRNA with single walled 
CNT IT 1 311 . Using state of the art all atom molecular dynamics simulation along with the ab-initio 
quantum mechanical calculations, we give a comprehensive understanding of the structure and 
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thermodynamics of the siRNA-graphene complex. In the supplementary material |37[|, we dis- 
cuss the effect of force fields (FF) on the unzipping and adsorption of siRNA/dsDNA on graphene. 



II. COMPUTATIONAL DETAILS 



A. Classical simulations 



We have used AMBER1 1 suite of programs 113 811 for simulating the systems with Amber 2003 
(along with ff99) force fields Ii39il and the TIP3P model I140n for water. Latest improvements 
of torsion angle parameters are reported for DNA/RNA in a new force field called parmbscO 
and fflO I14ll44-3ll . The comparison of force fields ff99, parmbscO and ff 10 are discussed in the 
supplementary material. Similar conformational changes are observed with ff99, parmbscO and 
ff 10. Use of ff99 also helps us to compare the present results with our previous simulations on 
siRNA-dendrimer complex that used ff99 Il29ll . The initial structure of the small interfering 
ribonucleic acid (siRNA) was taken from the protein data bank (PDB code: 2F8S) ll44n . The 
sequence of the siRNA used is r(UU AGA CAG CAU AUA UGC UGU CU) 2 with sticky ends of 
sequence UU on the two ends of the strands. We have built graphene sheet of 140 x 140 A 2 wider 
dimensions to ensure sufficient sliding area for the siRNA before optimum binding on graphene. 
For comparison, we have also simulated dsDNA of the same length on graphene. The dsDNA 
has same sequence as the siRNA where the nucleobase Uracil (U) is replaced by Thymine (T) 
nucleobase, i.e., d(TT AGA CAG CAT ATA TGC TGT CT) 2 . Structure of the dsDNA in B-form 
was generated using NUCGEN module of AMBER [38] . The siRNA/dsDNA-graphene complex 
structure is then solvated with TIP3P model water box using the LEaP module in AMBER 11 
113 811 . The box dimensions were chosen such that there is at least 20 A solvation shell from the 
surface of siRNA/dsDNA-graphene complex to the edge of the water box. In addition, some water 
residues were replaced by 44 Na + counterions for siRNA and 42 Na + counterions for dsDNA to 
neutralize the negative charge on the phosphate backbone groups of the siRNA/dsDNA structure. 
This has resulted in box dimensions of 188 x 186x72 A 3 with total system size of 214602 atoms 
for siRNA-graphene complex and 188x186x69 A 3 with total system size of 206865 atoms 
for dsDNA-graphene complex. The crystal structure of siRNA and the initial system containing 



siRNA and graphene with added water plus neutralizing counterions are shown in Figures 1(a) 
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and |l(b)[ respectively. 



We have modeled the carbon atoms in graphene as uncharged Lennard- Jones particles with LJ 
parameters listed in Table HI The bonded interactions viz., stretching, torsion and dihedral terms 
were also included. To keep the graphene fixed during simulations, all the carbon atom positions 
in graphene were restrained with harmonic potential of spring constant of 1000 kcal/mol-A 2 . The 
translational center of mass motions were removed every 1000 steps. The system is subjected to 
standard simulation protocol described in Refs. HQ- 

The long range electrostatic interactions 



were calculated with the Particle Mesh Ewald (PME) method Il47h using a cubic B-spline 
interpolation of order 4 and a 10~ 5 tolerance is set for the direct space sum cutoff. A real 
space cutoff of 9 A was used both for the long range electrostatic and short range van der Waals 
interactions with a non-bond list update frequency of 10. The trajectory was saved at a frequency 
of 2 ps for the entire simulation scale of 55-85 ns for each system. We have used periodic 
boundary conditions in all three directions during the simulation. Bond lengths involving bonds 



to hydrogen atoms were constrained using SHAKE algorithm [|48|] . This constraint enabled us 
to use a time step of 2 fs for obtaining long trajectories of each 50 ns. During the minimization, 
the siRNA/dsDNA-graphene complex structures were fixed in their starting conformations using 
harmonic constraints with a force constant of 500 kcal/mol-A 2 . This allowed the water molecules 
to reorganize which eliminates bad contacts with the siRNA/dsDNA and the graphene. The 
minimized structures were then subjected to 40 ps of MD, using 1 fs time step for integration. 
During the MD, the system was gradually heated from to 300 K using weak 20 kcal/mol-A 2 
harmonic constraints on the solute to its starting structure. This allows slow relaxation of the 
siRNA/dsDNA-graphene complex structure. Subsequently, simulations were performed under 
constant pressure-constant temperature conditions (NPT), with temperature regulation achieved 
using the Berendsen weak coupling method \a9\ (0.5 ps time constant for heat bath coupling and 
0.5 ps pressure relaxation time). Constant temperature-pressure MD was used to get the correct 
solvent density corresponding to experimental condition. Finally, for analysis of structures and 
properties, we have carried out 50 ns of NVT production MD with 2 fs integration time step using 
a heat bath coupling time constant of 1 ps. 

The binding free energy for the non-covalent association of two molecules in solution can be 
written as AG (A +B — > AB) = Gab — Ga~ Gb- For any species on the right hand side G(X) = 
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H(X) - TS(X), accordingly, 

AGbmd = bHbind ~ TASbj n d (1) 

where AZ/bind = AEgas +AG so i. Here Af/bind is the change in enthalpy and is calculated by summing 
the gas-phase energies (A£g as ) and solvation free energies (AG so i); E gas = E e \ e + £vdw + E[ nU 
where i?eie is the electrostatic energy calculated from the Coulomb potential, £Vdw is the non- 
bond van der Waals energy and E- m is the internal energy contribution from bonds, angles and 
torsions. G so i = G es + G ne s where G es is the electrostatic energy calculated from a Generalized 
Born (GB) method and G nes is the non-electrostatic energy calculated as y x SASA + 13; where y 
is the surface tension parameter (y= 0.0072 kcal/mol-A 2 ), SASA is the solvent-accessible surface 
area of the molecule and [3 is the solvation free energy for a point solute_((3 = 0). For the entropy 
calculation we have used two-phase thermodynamic (2PT) model [5(1 52, 53] , based on density 



of states (DoS) function. The DoS function can be calculated from the Fourier transform of the 
velocity auto-correlation function which provides information on the normal mode distribution of 



the system. The method has found successful application in several related problems [50-54]. 



6. Quantum simulations 

We have carried out quantum chemical calculations to understand the interaction of the planar 
nanographene with the nucleosides of RNA and DNA. For the present investigation, we have 
taken a graphene sheet of (6 x 6) dimension with alternate armchair and zig-zag edges with C-C 
bond lengths of 1.42 A, C-H bond length 1.09 A and all the angles kept at 120°. The edges are 
terminated with hydrogen atoms to avoid any unwanted terminal effects [55] • We have modeled 
adenine, guanine, cytosine, thymine and uracil with furanose sugar connected to the respective 
nucleobases by the [3-glycosidic bond. As structural features that distinguish siRNA from DNA 
are the presence of uracil nucleobase and hydroxyl-OH group in the constituted ribose sugar in 
siRNA, we have modeled thymine with the deoxyribose sugar and all other bases are modeled 
with ribose sugar to be consistent with siRNA structure. These structures act as miniature models 
of siRNA and DNA to understand the interaction with the planar graphene molecule. Initially the 
nucleosides are placed vertically above, around 4 A from the center of the planar hexagonal ring 
of the graphene and the nucleobases lie parallel to the graphene sheet. All the model building was 



done with the help of Discovery studio 2.0 Q56Q and Molden [57] software 
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We have optimized all the graphene-nucleosides complex systems to get the energy 
minimized configuration. For all the quantum chemical calculations, we have employed dis- 
persion corrected density functional approach using (OB97XD/6-31G** basis setQ58|] using 
Gaussian 09 Ii59il . The total interaction energy of each system has been calculated using 
E int = E (complex) — EX (isolated graphene) — EY {isolate d nucleoside)-!- BSSE, where BSSE 
represents basis set superposition error, arising from the overlapping of the atomic orbitals. The 
BSSE corrections has been calculated using Boys and Bernardi function counterpoise method 
J^oll . We have also carried out frequency calculation on the optimized geometry using the same 
method and basis set. To analyze the degree of binding, we have also carried out the charge 
transfer analysis of the nucleosides and graphene sheet using the Natural Bond Orbital (NBO) 



approach [61 



63] 



III. RESULTS AND DISCUSSION 



A. siRNA Unzipping 



Figure [T(c)1 shows the instantaneous snapshots of siRNA as it binds to graphene. The unzipping 
of base-pairs starts within 3 ns and by 22 ns, 41 of the initial 48 Watson-Crick (WC) H-bonds are 
broken resulting in almost complete unzipping of siRNA into two single strand RNAs. Figure [2] 
shows the number of intact H-bonds of siRNA and dsDNA with time when bound to graphene. 
Breaking of WC H-bonds in siRNA/dsDNA manifest the deformation mechanism of nucleic acid 



molecule 1164146611 . It can be seen that a large decrease of H-bonds occurs in the first few nanosec- 
onds when siRNA starts unzipping and within 22 ns, 41 H-bonds are broken leaving only 7 intact 
H-bonds. The unzipped bases are then free to interact with graphene surface via van der Waals 
forces. Figure [2] shows that H-bonds are saturated at 7. At this stage, siRNA optimally binds to 
graphene and the complex is very stable after 22 ns.In the optimum bound configuration where 
H-bonds are constant with time, there are few transient H-bonds due to room temperature ther- 
mal fluctuations. Because of these transient H-bonds, the curve has fluctuations about its mean 
value with standard deviations ranging from 1 to 3 H-bonds. The different nucleosides interact 
with different strength with graphene making the binding possible. Interestingly, dsDNA of same 
sequence and length show much less unzipping and binding on graphene compared to siRNA. 
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Within 6 ns, only 8 of total 48 H-bonds in dsDNA get unzipped and remain constant at 40 H- 
bonds throughout the rest of simulation time of 44 ns. This difference has its origin in the extra 
hydroxyl group of uridine that is present in siRNA. To have a molecular level understanding of 
this binding affinity we calculate the binding free energy of different nucleosides with graphene 
both from the classical MD simulations (using MM-GBSA method) as well as from our dispersion 
corrected DFT calculations. 

B. Binding free energy 

Figure |3] shows the enthalpy contribution to the total binding free energy as a function of time 
as siRNA binds to graphene. In the plot, we have marked the time interval at which optimal 
binding happens as reflected by the binding energy. After this binding the complex is stable for 
the entire duration of the simulation with fluctuations ranging only 1.5 % of its average value in 
the optimum bound state. From the stable trajectory of siRNA-graphene complex, the enthalpy 
contribution (A/^md) to the total binding free energy was calculated for 250 snapshots separated 
each by 2 ps. The enthalpy contribution (AHbi lu i) to the total binding free energy is -562.6 ± 
6.2 kcal/mol. Entropy is calculated at every 5 ns along the trajectory. For this, we simulate 
the system for 40 ps with velocities and coordinates saved at a frequency of 4 fs. The velocity 
auto-correlation function converges within 10 ps. When siRNA is binding with graphene during 
initial stage, siRNA entropy decreases since the graphene substrate suppress the fluctuations of 
unzipped bases (inset of Figure [3]). After siRNA binds optimally to graphene, the entropy starts 
increasing again due to inherent kinetics of unzipped bases. However the entropy contribution 
(TASbind) to the total binding free energy is small compared to the enthalpy contribution (AHbind) 
arising due to dispersive interaction between the graphene and the aromatic nucleobases. The 
enthalpy and entropy contributions in Eqn. [Ogive the total binding free energy (AGbind) of siRNA 
when binding to graphene. The value of AGu n d is -573.0 ± 8 kcal/mol. However, the binding 
energy of the dsDNA bound to graphene is calculated to be -190 ± 9 kcal/mol which is very low 
compared to siRNA bound to graphene. 

To get further molecular level picture of the binding mechanism we also compute the 
histograms of the closest approach of nucleoside to graphene. This will allow us to understand 
the relative binding affinity as well as how different nucleosides are oriented on the graphene 
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substrate. We track the center of mass position of different nucleobases with respect to graphene 
as a function of simulation time. We calculate the closest approach of nucleoside to graphene 
as the perpendicular distance r± = z • AB^ = B z — A z (because graphene is in x — y plane) 
where A and B are the centers of mass position of nucleoside and graphene, respectively. r± 
is calculated for 50 ns trajectory and the histogram of nucleosides is shown in Figure |4(a)| 
From this plot, we can understand the probability (P(r±)) of finding the nucleoside at a given 
position r± from the graphene. Using histogram, the free energy of nucleoside is calculated as 
F(r±) = —kgT ln(P(r±)) (where kg is the Boltzmann constant) and plotted as a function of 
distance in Figure |4(b)[ The minimum in the free energy indicates the most optimum bound 
configuration for nucleoside-graphene complex. In the inset we show the zoomed part of the 
occurrence of free energy minima. F(r±) has minima at 3.775 A , 3.790 A , 3.795 A , 3.840 
A and 3.860 A for guanosine, Thymidine, adenosine, uridine and cytidine, respectively. Guanine 
can form strongest H-bonding interaction and stacking interaction with graphene. The distance of 
closest approach in the increasing order for guanosine, thymidine, adenosine, uridine and cytidine 
is r op t(G) < r op t(r) ~ r op t(A) < r opt (U) < r opt (C). Hence guanine has most interaction strength 
with graphene and cytidine has least interaction strength asG>T~A>U>C. The binding 



14] 



energy order of nucleosides is consistent with experimental and theoretical calculations 1112 - 
The stable complex structure with most of the base-pairs already unzipped in siRNA can be 
delivered to the target virus infected cell for RNAi applications. Since siRNA has to undergo 
unwinding process with the effect of RISC, our proposed delivery mechanism by graphene 
possesses potential advantages in achieving RNAi. Toxic effects of graphene inside cell may be 



suppressed with proper surface functionalization Q67 



69|]. 



C. Structural deformation 

Snapshots shown in Figure Q] indicate that the siRNA molecule exhibit large structural 
deformation on binding to graphene. This structural deformation is characterized by the number 
of siRNA atoms that come close to the graphene in a specified cutoff distance and the root 
mean square deviation (RMSD) of siRNA with respect to its crystal structure. We calculate 
close contacts N c when any siRNA atoms are within 5 A of the graphene sheet. The number 



of close contacts N c between siRNA and graphene is plotted in Figure |5(a)| as a function of 
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time. Since siRNA is getting unzipped within 21 ns, N c is increasing rapidly within 21 ns and 
reaches a constant value of 710 siRNA atoms. The fluctuations are very less in N c after the 
complete binding of siRNA to the graphene. Interestingly, for dsDNA N c is much less than 
that of siRNA. The dsDNA has only 240 atoms within 5 A from graphene sheet in the stable 
configuration. This also demonstrates very low binding affinity of dsDNA compared to siRNA 
with graphene. In Figure |5(b)| we plot the RMSD of siRNA/dsDNA as a function of time for 
both the siRNA/dsDNA-graphene complex. For the calculation of RMSD, the reference structure 
of siRNA is taken to be the crystal structure of siRNA after initial minimization. Note that the 
plot shows RMSD for only production NVT simulation time scale. In the most optimum bound 
configuration, the average RMSD ((RMSD)i, ouru j) of siRNA is 19.4 ± 0.4 A on graphene whereas 
(RMSD) bound of dsDNA is 8.5 ± 0.6 A. As the binding of siRNA is more on graphene, the siRNA 
structure deforms leading to a large value of (RMSD)b ounc {. 



D. Unzipping kinetics 

In Figures [6(a)] and [6(b)] , we plot the number of contacts for siRNA as well the distance between 
the two strands of the double stranded siRNA (ds- separation) at three different temperatures. It 
takes few ns to nucleate the unzipping. Once critical numbers of contacts are created between the 
siRNA and graphene, unzipping starts from one end (end2 for the current situation). Unzipped 
bases at one end help to make more contacts with graphene and thereby enhancing the interaction 



between siRNA and graphene. This facilitates unzipping at the other end. Figure [6(a)] shows that 
it takes almost 3 ns at 300 K and less than 1 ns at 340 K for the critical number of contacts to 
be created. Once this is done, strong interaction between graphene and siRNA force the rapid 
unzipping of siRNA as is evident from the rapid increase of ds-separation in Figure |6(b)[ To 
get an estimate of the unzipping time (x), we plot the unzipping probability distributions (fhb) at 
three different temperatures in Figure [7] and fit them to single exponential functions as done in 
ref. [70]. This gives rise to x = 9.8 ns, 8.4 ns and 5.0 ns at temperatures 290 K, 300 K and 340 
K, respectively. By fitting the unzipping time as a function of temperature to x = %oe Ea ^ kBT , we 
get Xo to be 100.8 ps and the activation energy, E a to be 2.637 kcal/mol or 0.114 eV. Plot of Inx 
versus ^ and the fit was shown in the inset of Figure [7] 
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E. Insights from Quantum simulations 



To understand the difference in binding affinity of siRNA and dsDNA with graphene we have 
calculated the binding energy of the graphene/siRNA and graphene/dsDNA miniature complexes 
using dispersion corrected DFT method (DFT-D). The sequence of the siRNA studied has A:U 
base pairs at both the ends, which is also the case of most of the confirmed siRNA sequences [ 7lj . 
which are known to open up quite easily as compared to G:C base pairs. It may be noted that most 
of the DNA/RNA oligonucleotide sequences whose three-dimensional structures in double helical 
forms are available have G:C base pairs at both the termini Ii72ll . The siRNA sequences need 
to unzip soon for their functionality and probably that drives design of sequences with terminal 
loosely bound A:U base pairs. However, it appeared that the dsDNA molecule, containing 
terminal A:T base pairs does not unzip at physiological condition. In order to check whether this 
is due to stronger attraction in A:T as compared to A:U base pairs, we have optimized both these 
base pairs using wB97XD/6-31G** method and found their interaction energies to be -15.80 and 
-15.91 kcal/mol, respectively. This clearly indicates that the terminal A:U base pair is not weaker 
one as compared to its DNA counterpart. Thus, the other component of differential interaction, i.e. 
interactions between Uracil residues and Thymine residues with graphene might be driving the 
RNA molecules to unzip. We are mainly interested in thymidine-graphene and uridine-graphene 
complex systems since these are the principal nucleosides, which can differentiate between DNA 
and RNA. Moreover these nucleosides remain unpaired in the siRNA as well as dsDNA. On 
analyzing the optimized geometry of the complex systems as shown in Figure [8] (graphene-uridine 
and graphene-thymidine complex), we find that in both the cases 03'-H3' of the constituent sugar 
points towards the planar graphene sheet with close approach forming 0-H...7T contacts. In case 
of thymidine nucleoside, the closest OS'-rD'^-ring center is found to be around 2.54 A , and the 
angle <03'-H3'...ring center> is obtained as 129.35°. These distances and angles are sufficient 
to form O-H...71; types of H-bonds. While in case of uridine nucleoside, the 03 / -H3 / ...ring center 
bond distances are found to be 2.34 A , and the angle <03 / -H3 / ...C> is obtained as 165.32°, 
forming significantly stronger H-bond between the graphene and uridine sugar, as compared 
to that of thymidine sugar complex \T^\ . The 02 and 04 groups of the thymine and uracil 
also interact with the graphene sheet, but the magnitudes of interaction seems to be very low as 
compared to that of 0-H...7C types of contact. Nevertheless these carbon oxygen atoms can form 
lone pair. . . 7t type of contacts giving extra stabilization I174n . In addition to these contacts, the 02' 
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may also form H-bond with the graphene, which however was not found in the energy minimized 
structures, but can not be ignored at physiological temperature and between the oligonucleotides. 
Furthermore several of 2'-OH groups, which are equivalent to 3'-OH, are present in the siRNA 
strands and absent in the dsDNA strands. Thus the siRNA strands tend to dissociate from their 
double helical forms and bind strongly to the graphene sheet. 



The BSSE corrected interaction energy has been calculated for all the complex systems 
and they follows the trends G>A>U>T>C. Interaction energy of the graphene-uridine 
nucleoside is found to be -22.26 kcal/mol, while that of graphene-thymidine nucleoside is found 
to be -20.30 kcal/mol. So we can infer that uridine nucleoside interacts with the graphene more 
strongly than that of the thymidine nucleoside. These interaction energy values also well correlate 
the H-bond lengths and angles values obtained in both the cases. Our previous studies Ii36ll on 
interactions of the graphene with the nucleobases gives the interaction energy strengths as G > A 
> C > T > U. So it proves that inclusion of sugar in the nucleobases alters the interaction energy 
strengths, therefore plays significant role in stabilizing the systems due to availability of more 
hydrogen bond donor sites. 



Frequency calculation of the entire complex as well as the isolated systems give no imaginary 
frequency indicating the structures are at their local minima. Frequency calculation enables 
us to carry out thermochemical analysis of the system. All the calculations are carried out at 
298.15 K and 1 atm. pressure. We have calculated the change in enthalpy and free energy of 
the systems. The AH of graphene-uridine is found to be -26.39 kcal/mol, whereas that of AH of 
graphene-thymidine is calculated to be -23.80 kcal/mol. Similarly the AG of graphene-uridine and 
graphene-thymidine are found to be -1 1.95 and -10.22 kcal/mol, respectively. As it is well known 
that the system is more favorable with increase in the negative value of AG, graphene-uridine 
complex is more stable. The formation of stable graphene-uridine nucleoside complex may 
initiate and enhance the unzipping of the siRNA structure as observed by our counterpart MD 
simulation studies. 



We have also calculated the NBO charges of the thymidine and uridine nucleosides of the 
graphene-nucleosides systems and compared with the isolated nucleosides. The difference in 
NBO charges of the major hydrogen bond donor atoms of the nucleosides which interacts with 
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the graphene sheet are given in Table Hfl We observed that charge transfer is more significant 
for the 03'-H3' and 04 atoms for uridine molecule with graphene, since they interacts strongly 
with the graphene sheet. The 02 of thymidine shows negligible amount of charge transfer with 
the graphene sheet. It may also be noted that due to the methyl group of thymine, close to the 
04 atom, the 04 atom of thymine may not be allowed to come close to the graphene plane, 
thereby reducing its interaction strength. We can conclude that uridine interacts with the graphene 
sheet more strongly than thymidine, which is well correlated with the hydrogen bond strengths, 
interaction energy, and thermochemical analysis of the systems. 



IV. CONCLUSION 



We demonstrated very unusual phenomena of complete siRNA unzipping and binding on 
graphene substrate. One of the major goals of the current study is to understand the binding 
mechanism of siRNA/dsDNA on graphene. Our study also shows that siRNA unzips and binds to 
graphene forming graphene-siRNA hybrid. This allows us to study the siRNA unzipping kinetics. 
siRNA unzipping is very important in the context of RNAi therapeutics where a short siRNA 
enters into cell and gets unzipped for its further action. Studying unzipping kinetics through 
graphene also offers an alternate route to nanopore assisted unzipping where an electric field is 
applied to translocate dsRNA/dsDNA. The stable graphene-siRNA hybrid may also be used for 
efficient delivery of siRNA. The complex may penetrate the hydrophobic regions of the bilayer 
due to hydrophobicity of graphene. Inside the cell, the graphene-siRNA may remain as complex 
between graphene and two single stranded RNA chains. One of the chains might easily dissociate 
from the complex whenever a competitive messenger RNA chain approaches the complex. This 
can silence the required gene. Another interesting observation is that dsDNA of same sequence as 
siRNA except thymine in place of uracil has less unzipping and less binding on graphene. This 
interesting property could be used to detect or separate siRNA and dsDNA molecules. We support 
these findings through long classical MD simulation as well as calculating the relative binding 
affinity of nucleosides with graphene through dispersion corrected DFT methods. It is shown 
that the unpaired uracil residues make strongest contacts with the graphene molecule through 
van der Waals and specific H-bonding interaction involving 2 / -0H group of the ribose sugar. 
These interactions can be responsible for the double helical siRNA to unzip, which are stabilized 
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subsequently by several such 0-H...71 interaction. The equivalent double stranded DNA does not 
have the -OH group and hence remains stable throughout the simulation time. The spontaneous 
unzipping helps us to study the siRNA unzipping kinetics for the first time. Unzipping time is of 
the order of 5-10 ns and decreases with increasing temperature. Unzipping time follows Arrhenius 
behavior and allows us to get an estimate of the energy barriers for the siRNA unzipping. In 
contrast to the unzipping kinetics study through nanopore unzipping which requires application 
of voltage and takes longer time 117 511 . unzipping in graphene is very fast and happen spontaneously. 
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Figures 




FIG. 1: (a) siRNA crystal structure (pdb code 1F8S) and (b) the initial simulation system setup where 
siRNA-graphene complex was solvated with water and neutralizing Na + counterions. (c) Snapshot of 
siRNA on graphene after optimum binding, (d) and (e) shows the zoomed part of the completely un- 
zipped siRNA strands whose bases are in ji — n interaction with graphene. The water and counterions were 
not shown in (c), (d) and (e) for image clarity. The snapshots were rendered using VMD software package 
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FIG. 2: Number of intact Watson-Crick H-bonds in siRNA and dsDNA as a function of time with graphene. 
The siRNA has lesser intact WC H-bonds than dsDNA. 
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FIG. 3: Enthalpy of siRNA when binding to graphene as a function of time. Correspondingly, the entropy 
of siRNA is shown in the inset. The enthalpy gets saturated within 18 ns to -562.6 kcal/mol. In the most 
optimum bound configuration, the binding free energy that includes enthalpy and entropy contributions is 
-573.0 ± 8.0 kcal/mol. 
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(a) 



(b) 



FIG. 4: Nucleosides interaction with graphene from classical MD: (a) Histogram of nucleosides position 
from graphene surface. In the plot r± is the perpendicular distance between the center of masses of nu- 
cleoside and graphene as shown, (b) Free energy of the nucleoside-graphene complex as a function of 
r±. Minima in the free energy curve indicates the most optimum bound configuration for the nucleoside- 
graphene complex. Inset is the zoomed part around free energy minima. The minima occurs at a value in 
the increasing order for guanosine, thymidine, adenosine, uridine and cytidine. The smaller values of r± 
may imply a larger interaction strength with graphene. Hence guanine has most interaction strength with 
graphene and in decreasing order thymidine, adenosine, uridine and cytidine. 
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FIG. 5: Number of the close contacts N c and RMSD of siRNA and dsDNA as a function of time. For 
siRNA, N c and RMSD are maximum compared to dsDNA. 
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FIG. 6: Temperature effects: (a) Number of close contacts of siRNA on graphene (b) ds-separation of 
siRNA while unzipping from both the ends. 
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FIG. 7: Unzipping probability denoted by the fraction of broken WC H-bonds, as a function of time for 
different temperatures. 
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(a) 




(b) 

FIG. 8: Quantum DFT-D Optimized geometry of the complex systems: (a) graphene-uridine nucleoside and 
(b) graphene-thymidine nucleoside. 



Tables 



TABLE I: Lennard- Jones interaction parameters used for carbon atoms in graphene and counterions Na + . 



Atom 


e/k B (K) 


a (A) 


C 

Na+ 


43.7 
1.41 


3.40 
3.328 



TABLE II: Non bonded orbital (NBO) charges of the thymidine and uridine nucleosides 



Atom No. 


Graphene + Thymidine 


Graphene + Uridine 




difference in NBO charge 


difference in NBO charge 


03'-H3' 


0.005 


0.007 


02 


0.002 


0.000 
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04 



0.007 



0.014 



23 



43: 
a* 



Of 



Supplementary material for 
"Unraveling siRNA Unzipping Kinetics with Graphene" 

Santosh Mogurampelly 1 , Swati Panigrahi 2 , Dhananjay 
Bhattacharyya 2 , A. K. Sood 3 , and Prabal K. Maitir 

1 Centre for Condensed Matter Theory, Department of Physics, 
Indian Institute of Science, Bangalore 560012, India 
2 Biophysics Division, Saha Institute of Nuclear Physics, Kolkata 700064, India 
3 Department of Physics, Indian Institute of Science, Bangalore 560 012, India 



*To whom correspondence should be addressed; Electronic address: 



1 



maiti@physics.iisc.ernet.in 



I. FORCE FIELD (FF) EFFECTS 



The results of siRNA unzipping mechanism on graphene discussed in the main paper are based 
on ff99 force field yj] . The ff99 force field has been widely used in the last decade to study the 
structure and thermodynamics of nucleic acids. However, ff99 is found to produce unphysical 
torsion angle transitions in 'oc/y' in longer MD simulations. These irreversible 'oc/y' transitions 
introduce severe distorsions in nucleic acids. Based on these torsion angle parameterization 
artefacts, we want to answer if the large structural transitions in siRNA/dsDNA when adsorbed to 
graphene are due to torsion angle artefacts of ff99. Recently, torsion angle parameters for 'oc/v' 
and glycosidic dihedral of nucleic acids have been refined in parmbscO and fflO [4] 

force fields, respectively. Note that ff 10 is equivalent to ff99 + parmbscO for DNA and ff99 + 
parmbscO + % for RNA. 



We have studied adsorption of siRNA/dsDNA on graphene with ff99 and parmbscO. In 
addition, siRNA adsorption on graphene has been studied with ff 10 also and as mentioned for 
dsDNA, parmbscO is equivalent to using fflO. Snapshots and plots of siRNA/dsDNA on graphene 
with various force fields were shown in Figures \T\ and [21 respectively. The major conclusion 
from these studies is that with parmbscO and fflO also, unzipping and strong adsorption of 
siRNA on graphene is observed. In comparison, dsDNA has less unzipping and adsorption 
on graphene compared to the unzipping and adsorption of siRNA. Some of the earlier studies 
have also demonstrated that ff99 parameters are as good as parmbscO set in describing the 
conformational richness of RNA |5kZ|] . In our simulation with fflO, less unzipping of siRNA is 
observed compared to ff99/parmbsc0 (Figure [2(a)] ) . In contrast, with parmbscO, dsDNA has more 



unzipping and adsorption than ff99 as it is clear from snapshots (Figure Q]) and WC H-bonds, 
close contacts and RMSD (Figure [2]). In spite of these differences, the qualitative features of the 
siRNA unzipping and adsorption on graphene have been reproduced with all the force fields under 
consideration. Hence, we conclude that the large structural changes observed with ff99 force field 
are not because of the artefacts of torsion angle parameters. 
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FIG. 1: FF effects: Snapshots of siRNA and dsDNA unzipping and adsorption on graphene using various 
force fields. The dimension of graphene sheet is 140x140 A 2 but only the portion of siRNA/dsDNA 
adsorption is shown in the figure. These snapshots are taken in the most optimum bound configuration 
both for siRNA and dsDNA. With ff 10 the unzipping of siRNA as can be quantified from the number of 
intact WC ITbonds (Figure |2(a)| ) is less but unzipping and adsoiption of siRNA on graphene is qualitatively 
similar to that observed with ff99/parmbsc0. However, parmbscO produces more unzipping of dsDNA than 
££99 force field (Figure [2(a)]) . 
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FIG. 2: FF effects: Time series of the number of intact Watson-Crick H-bonds, close contacts and RMSD 
of siRNA/dsDNA while adsorbing on graphene for various force fields. These structural quantities verify 
that the unzipping and adsorption of siRNA/dsDNA on graphene can also be observed using parmbscO and 
ff 10. It is interesting to note that RMSD is same for both the siRNA and dsDNA with parmbscO while they 
are very different when using ff99/ffl0. 
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